Haciendo Mapas con Python

Stephanie Orellana

Pycon Chile 2022

Sobre mi

> Stephanie Orellana (@sporella/@sporella.dev)

> Ingeniera Agrónoma y Magister en Recursos Naturales

> Científica de Datos

> PyLadies Valparaíso

Mapas :)

Tipos de datos espaciales

Vectoriales

  • .shp (shapefile)
  • .geojson
  • .csv
  • …

Raster

  • .tif
  • .nc (netcdf)
  • .asc (ascii)
  • .hdf
  • …

Tipos de datos espaciales (ejemplos)

Vectoriales

  • Límites administrativos
  • Coordenadas puntuales
  • Caminos
  • Tracks de GPS

Raster

  • Imágenes satelitales
  • Interpolaciones espaciales
  • Modelos climáticos
  • Modelos de elevación digital
  • Imágenes de drones

Paquetes

Vectoriales

> geopandas

Raster

> xarray

> rasterio

Visualización

> matplotlib

> folium

> plotly

> bokeh

> cartopy

Pero, ¿dónde consigo los datos?

  • IDE (Infraestructura de datos espaciales)

  • Webs gubernamentales

  • NASA / COPERNICUS

  • Google Earth Engine

  • Open Street Map

  • Datos de campo (GPS)

Fuentes de Información

En Chile

Mundo

geopandas

GeoPandas, amplia la librería de ciencia de datos pandas añadiendo soporte para datos geoespaciales.

Nuestro primer mapa

import geopandas as gpd
comunas = gpd.read_file("data/comunas_valparaiso_geo.shp")
comunas.plot("Provincia", legend = True, figsize = (6, 6));

Personalizándolo con matplotlib

Código 
import matplotlib.pyplot as plt
from matplotlib_scalebar.scalebar import ScaleBar
import geopandas as gpd

comunas = gpd.read_file("data/comunas_valparaiso_geo.shp")

d = comunas.to_crs(32719)

ax = d.plot("Provincia", legend = True, cmap = "Pastel2", 
legend_kwds={'loc': 'upper left', 'bbox_to_anchor':(1, 0.5)}, 
figsize = (6, 6))


ax.set_axis_off()

ax.add_artist(ScaleBar(1))

x, y, arrow_length = 0.1, 1, 0.1
ax.annotate('N', xy=(x, y), xytext=(x, y-arrow_length),
            arrowprops=dict(facecolor='black', width=2, headwidth=8),
            ha='center', va='center', fontsize=10,
            xycoords=ax.transAxes);

geopandas + folium

 

Tal como el método plot(), geopandas tiene un método llamado explore() que nos permite hacer visualizaciones dinámicas con folium.

folium nos permite utilizar en Python la librería leaflet.js

 

import geopandas as gpd

paradas = gpd.read_file("data/paradas_transantiago.json")
paradas.explore()

geopandas .explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

folium plugins

from folium import plugins
from folium.plugins import HeatMap

# guardamos nuestro mapa inicial en un objeto
m = paradas.explore()

# nos aseguramos que las coordenadas presentes en la tabla sean de tipo numérico
paradas['stop_lat'] = paradas['stop_lat'].astype(float)
paradas['stop_lon'] = paradas['stop_lon'].astype(float)

# generamos una lista con latitudes y longitudes
heat_data = [[row['stop_lat'], row['stop_lon']] for index, row in paradas.iterrows()]

# Aplicamos la función HeatMap e incluimos al mapa inicial
HeatMap(heat_data).add_to(m)
m

Densidad de puntos

Make this Notebook Trusted to load map: File -> Trust Notebook

Bajar datos desde OpenStreetMap

Podemos descargar datos desde Open Street Map con el paquete osmnx

import osmnx as ox

lugar = "Renca, Chile"
tags = {'landuse': True}   

uso_suelo = ox.geometries_from_place(lugar, tags)
uso_suelo.head()

Bajar datos desde OpenStreetMap

source geometry name nodes landuse operator plant:method plant:output:electricity plant:source power start_date addr:housenumber addr:street denomination ways type man_made
element_type osmid
way 25291350 NaN POLYGON ((-70.73081 -33.40046, -70.73064 -33.4... NaN [275503907, 275503908, 275503909, 275503907] village_green NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
26409096 https://generadora.cl/ POLYGON ((-70.69016 -33.41707, -70.69004 -33.4... Central Termoeléctrica Renca [607764617, 621925434, 289338330, 289338012, 2... industrial Generadora Metropolitana combustion 479 MW gas;oil plant 1962 NaN NaN NaN NaN NaN NaN
26409384 NaN POLYGON ((-70.68475 -33.41905, -70.68359 -33.4... NaN [289342010, 289342012, 289342013, 289342014, 2... industrial NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
26409392 NaN POLYGON ((-70.68144 -33.41881, -70.67846 -33.4... Megacentro [289342102, 289342103, 289342104, 289342105, 2... industrial NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
26409395 NaN POLYGON ((-70.68224 -33.42129, -70.67755 -33.4... NaN [289342141, 289342142, 289342143, 289342144, 2... industrial NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

Y hacer nuestro mapa con explore

import geopandas as gpd
import osmnx as ox

lugar = "Renca, Chile"
tags = {'landuse': True}   

uso_suelo = ox.geometries_from_place(lugar, tags)

uso_suelo.explore(
  "landuse",
  tiles = "CartoDB Positron",
  tooltip = False,
  popup = ["landuse"],
  highlight = True
)

Y hacer nuestro mapa con explore

Make this Notebook Trusted to load map: File -> Trust Notebook

Rasters

Existen varias formas de visualizar rasters, pero acá mostraré la forma más simple que es usando matplotlib

import cartopy
import matplotlib.pyplot as plt
import xarray as xr
import rioxarray

r = xr.open_rasterio("data/soil_moisture.tif")

fig, ax = plt.subplots(dpi = 300, figsize=(4, 4), 
subplot_kw = dict(projection = cartopy.crs.Mercator())) 
r.plot(cbar_kwargs={'label': "Humedad de Suelo [mm]", 
                    "orientation": "horizontal"}, 
                    cmap = "PuBuGn")
plt.title("")
plt.axis('off')
plt.show()

Rasters

Palabras finales

  • ¡Hagan mapas!

  • Mapas son una representación en 2d de un mundo en 3d

  • Sesgos, ética

  • Calidad de los datos

Recursos

Libros

Cursos y tutoriales